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> Abstract 

o 

CN Rubber components and soft tissues are often subjected to large 

bending deformations "in service". The circumferential line elements 
• on the inner face of a bent block can contract up to a certain critical 

Q stretch ratio Acr (say) before bifurcation occurs and axial creases ap- 

pear. For several models used to describe rubber, it is found that 
Acr = 0.56, allowing for a 44% contraction. For models used to 
^ describe arteries it is found, somewhat surprisingly, that the strain- 

^ stiffening effect promotes instability. For example, the models used 

5-H for the artery of a 70 year old human predict that Acr = 0.73, allowing 

only for a 27% contraction. Tensile experiments conducted on pig skin 
indicate that bending instabilities should occur even earlier there. 

Keywords: large bending, nonlinear elasticity, bifurcation, strain-stiffening 
effect, soft tissue modeling. 
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1 Introduction 



When an elastic block is subjected to severe bending, it is expected that its 
inner curved face will eventually become unstable. This phenomenon is well 
captured by the theory of incremental nonlinear elasticity, which can account 
for the appearance of small-amplitude wrinkles, aligned with the axial direc- 
tion of the large bending. For a block made of rubber, Gent and Cho (1999) 
argue that this buckling should occur once circumferential line elements on 
the bent face are contracted up to the critical stretch (Acr, say) of surface in- 
stability for a semi-infinite body. In that latter case, the neo-Hookean model 
predicts (Biot, 1963) that Acr = 0.54, allowing hne elements to be contracted 
by 46%. Analyses looking into the bending instability (instead of the surface 
instability) of neo-Hookean blocks show that the finite size of a block (as 
opposed to the infinite size of a half-space) introduces dispersion (and thus 
links the size of a block to the number of wrinkles) but has little effect on the 
amount of possible contraction, bringing Acr from 0.54 to 0.56, irrespective 
of the dimensions of the block. 

Now the neo-Hookean model is often used to describe solids undergoing 
finite but moderate deformations. In fact, as recalled in Section [2] below, it 
encompasses the most general incompressible solid of third-order elasticity, as 
well as the Mooney-Rivlin model for rubber, when the deformation is a plane 
strain such as pure bending. For large strains however, the polymers chains in 
elastomers and the collagen fiber bundles in soft tissues align themselves with 
the direction of greatest stretch and their limiting extensibility is strongly felt. 
This phenomenon is the so-called strain- stiffening effect. In biological soft 
tissues the stiffening occurs at markedly lower strain levels than in rubber- 
like solids: broadly speaking, rubbers can be stretched at least 100%, whilst 
soft tissues can be stretched at most 100%, before they stiffen. 

There are two popular models of constitutive law to describe the strain- 
stiffening effect: the Gent model, which introduces a limiting chain parameter 
— and thus an upper bound for the range of possible stretches — , and the 
Fung model, which exhibits an exponential increase of stress with respect to 
strain — but no maximal stretch. The two models are presented in Section 
2, along with typical values for the stiffening parameters of soft tissues. In 
particular, we conducted tensile experiments on pig thoracic aortas and found 
good agreement with data published previously in the literature. 

Then in Sections 3 and 4, we show that when these two models are ad- 
justed to account for the actual physiological values of mammalian arteries. 
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they predict that bending instabihties appear early, at moderate amounts 
of bending. For example , the Gent model, with a stiffening parameter ad- 
justed to describe the thoracic artery of a 70 year old human, predicts that 
Acr = 0.73, which means that circumferential line elements can be contracted 
by 27% only, before instabihties arise. Further, the two models predict that 
stiffer tissues such as skin should become unstable almost as soon as they are 
bent. To estimate the stiffening parameters in those cases, we conducted ten- 
sile experiments on pig skin, and present the corresponding data and curve 
fitting results at the end of section (§4). 

These results suggest that other characteristics must be accounted for in 
the constitutive modeling of arteries and skin. In particular, other investi- 
gations e.g. (Destrade et al., 2008) show that the anisotropy displayed by 
biological soft tissues (as opposed to the isotropy of elastomers) plays a most 
significant role in stability studies. 

2 Basic equations 
2.1 Constitutive laws 

We are interested in hyperelastic strain energy densities {W = W{Xi, A2, A3) 
where the A's are the principal stretches) commonly used to model the be- 
havior of strain-stiffening incompressible solids. One of the simplest and best 
prototypes available to do this is the Gent model (Gent, 1996), 

2 \ Jm J 

where > is the shear modulus for infinitesimal strains and > is 
a stiffening parameter, associated with limiting chain extensibility (Horgan 
and Saccomandi, 2006). Hence a solid behaving according to this constitutive 
model cannot be stretched in uniaxial tension beyond a maximal stretch A^, 
say, which is the positive root of 

A^ + 2A-^- J^-3 = 0. (2) 

Typically, rubbers attain maximal strains ranging from 5 to 15; in contrast, 
biological soft tissues experience strain-stiffening effects much earher. For 
instance experimental data on young, healthy human (Horgan and Sacco- 
mandi, 2003) and bovine (Humphrey, 2003) arteries show that they can be 
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stretched up to a maximum of = 1-4, or even 2.0, while older, stiffer 
arteries can be stretched only up to 20% max (Am = 1.2). Accordingly, in 
what follows we take 20 < Jm < 200 as the typical range for rubbers, and we 
take 0.4 < Jm ^ 2.3 as the typical range for arteries. Other biological soft 
tissues belong to that latter range, including caterpillar muscle (Dorfmann 
et al., 2007), rabbit muscle (Davis et al., 2003), dura mater (Maikos et al., 
2008), brain tissue (Franceschini et al., 2006), etc. 

Another popular model in the biomechanics literature, which reflects a 
less pronounced strain-stiffening effect than the Gent model, is the Fung 
model, 

(3) 



2b 



gb(A^+Ai+Ai-3) _ ^ 



where > is the shear modulus at small strains and 6 > is the stiffening 
parameter. For a human young thoracic artery, b ~ 1.0, and for an older, 
stiffer artery, 6 ~ 5.5, typically (Horgan and Saccomandi, 2003). Accordingly, 
we take 1.0 < b < 5.5 as the typical range for arteries when using the Fung 
model. This range is also inclusive of Fung model fitting for human cerebral 
veins (6 ~ 1.7) and cerebral arteries (6 ~ 4.4), see Ho and Kleiven (2007). 

Using a Tinus Olsen tensile machine, we conducted tensile tests on several 
pig aortas, and found that their Gent and Fung parameters did indeed belong 
to those ranges, see Figure [T] for an illustrative example. 

Note that the Gent and the Fung models are both connected to the so- 
called neo-Hookean model, 

W = /i(A? + A^ + A^ - 3)/2, (4) 

in the limits Jm = oo, b = 0, respectively. 

In this connection we mention the so-called third-order elasticity model, 
often used to investigate the onset of nonlinear elasticity effects. When the 
general strain energy density of an incompressible solid is expanded in terms 
of powers of E, the Green strain tensor, it is found that at the third-order 
of truncation (Ogden, 1974; Hamilton et al., 2004), 

W = ijtY{E^) + lAtT{E^) , (5) 

where n is the infinitesimal shear modulus of second-order (linear) elasticity, 
and ^ is a Landau third-order elastic coefficient. Turning to a representation 
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stretch V stress (Fung model) 

0.16 r 




Figure 1: Tensile test on a pig thoracic aorta: nominal stress versus stretch 
ratio. The thick line corresponds to experimental data and the thin line, to 
the curve fitting using the Gent model (on the left) and the Fung model (on 
the right). Curve fitting analysis gives Jm — 1-4 with the Gent model and 
h ~ 1.3 with the Fung model. 



in terms of the principal stretches, it is found that at the same order, 

W={2^^+ ^) {XI + + - 3) - + ^) (A?A^ + A^A^ + A^A^ - 3) . 

(6) 

For a plane strain deformation such as pure bending, A3 = 1 at all times 
and then A2 = A|f ^ by incompressibility. In that case, comparison of (g and 
^ shows that the third-order elasticity model and the neo-Hookean model 
coincide. Note that in finite elasticity, ([6]) is the so-called Mooney-Rivlin 
strain-energy density, often used to describe rubbers. 



2.2 Pure bending of a block 

We start with a straight block of thickness 2A, width 2L, and height H, 
located in the region 

-A<Xi<A, -L<X2<L, 0<X3<H, (7) 

where Xi, X2, X3 are the coordinates in the rectangular coordinate system 
aligned with the edges of the block, giving the position of a material point 
in the reference configuration. 
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Figure 2: Pure bending of a thick block. Here the original length-to-thickness 
ratio is 2.0 and the bending angle is 7r/3. The circumferential line elements 
on the outer bent face are extended during bending; those on the inner face 
are contracted. Eventually, at a critical stretch of contraction, the inner face 
buckles, and axial wrinkles appear. 

By applying surface tractions on the the faces at X2 = ±L, we bend the 
block using the deformation (Green and Zerna, 1954) 

r = Vc? + 2X1/00, e = 00X2, z = X3, (8) 

where is a constant (to be determined later), a; is a prescribed constant, 
and r, 6, z are the cylindrical coordinates, giving the position of a material 
point in the deformed configuration. The bending angle is (p = 2ujL. The 
bent block is confined to the region 

Ta = \Jd - 2A/u <r <n = ^d + 2A/uj, -uL <e <ujL, < z < H, 

(9) 

where Va are Vb are the radii of the inner and outer faces. 
The corresponding principal stretches are 

Ai = {ior)-\ A2 = oor, A3 = 1, (10) 

showing that this is a plane strain deformation. These expressions also show 
that the strain energy density W{\i, A2, A3) is a function of only one variable. 
Choosing it as A = Ai, we call S the function defined by 

E{X)^W{X,X-\1). (11) 
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Then Rivlin (1949) shows that the principal stress component ai is given by 

ai = S(A) + ir, (12) 

where the constant K is to be determined from the boundary conditions. 

We assume that the bent faces at r = ra and r = Vb are free of normal 
tractions 

ai{ra) = (T,{n) = 0. (13) 

These two boundary conditions lead to the determination of the constants d 
in dSl) and in ( 12 ) for the strain energy densities (|l|), (g, and Q. 

Take for instance the neo-Hookean form (|4]): the boundary conditions 
read 

fi[icora)-^ + {coray-2] + 2K = 0, fx[{con)-^ + {unf -2]+2K = 0, (14) 
from which follows that (Green and Zerna, 1954) 



u ran = 1 

Then we find d and K as 

d = {l/u^)Vl+4u^A^, K = 
and the normal stress as 

,2 



(15) 



(^/2)[a;^(r^ + r,^)-2], (16) 



2r2 



(17) 



For the Gent and Fung models, d is also given by (16)i, whilst o"i(r) is equal 
to 



fiJr, 



In 



Jm ~t~ 2 



ran 



+ 2 - ^ - a 



and 



ii 
2b 



\rarf, J — p \ra J 



(18) 

respectively, see Kanner and Horgan (2008) and Demiray and Levinson (1982) 
for details. 

To complete the picture, we give expressions for a2, the principal Cauchy 
stress normal to the surfaces 6 = const., and for its moment A4 of the stresses 
on the faces 6 = ±uL about the origin: 



0"2 = d(rcri)/dr, M. = H ra2dr. 



(19) 



respectively. 
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3 Bending instability 



Instability in bending is triggered by the apparition of wrinkles/ creases on the 
inner face of the bent block. Their existence is governed by the incremental 
equations of equilibrium and of incompressibility (Ogden, 1984), 

div s = 0, div u = 0, (20) 

where s is the incremental nominal stress and u is the incremental mechanical 
displacement, and by the satisfaction of appropriate boundary conditions. 

In the present context of pure bending in nonlinear elasticity, these equa- 
tions have been written down several times (see (Coman and Destrade, 2008) 
and references therein) , usually as a second-order system of coupled differen- 
tial equations for the components of the displacement. Instead, we present 
them here as a first-order system for the displacement and the traction: this 
is the so-called Stroh formulation, which proves to be optimal for the sub- 
sequent numerical resolution of the boundary value problem. Omitting the 
details, we find that if the components u, v, w of the mechanical displacement 
and Srr, Sre, Srz of the traction are in the form 

{u, V, w, Srr, Sro, Sr.} = { [U{r), V {r) , 0, 5„(r), Sre{r),0] e^"^} , (21) 

where U, V, Srr, Sre, are complex functions of r only, then two equations in 



(20) are identically satisfied and there remains four independent equations. 



which can indeed be put in the Stroh form. The constant n in (21) is the 
circumferential number, and is determined from the condition that there are 
no incremental normal tractions on the end faces 6 = ±uL: this happens 
when (Haughton, 1999) 

n = p7i/{ujL), (22) 



for some integer p, which we call the mode number. Hence, the solution (21) 
describes the p creases appearing on the inner surface of the bent block. 

Introducing the four-component displacement-traction vector (Shuvalov, 
2003) 

rj = [U,V,irSrr,irSreY, (23) 
we find that the incremental equations of equilibrium can be arranged as 

-^r/(r) = -G{r)r}{r). (24) 
dr r 
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Here the matrix G has the following Stroh structure, 



G 



i 


—n 








n(l — ai/a) 


-i(l - (Ti/a) 





-1/a 


Kii 




— i 


—n{l — (Ti/a 


-i/«i2 




—n 


i(l - ai/a) 



(25) 



where 



Kii = 2(/3 + a — ai) + [7 — (a — o"i)^/a] 

K12 = n (2/3 + a + 7 - (t^/ a) , 

K22 = 7 - (« - o"i)^/q; + 2n^(/3 + a - ai), 



(26) 



and the quantities 7, a, and /3 are given in turn by (Dowaikh and Ogden, 
1990) 



7 = AS'(A)/(A^-1), 



a 



2(3 = A'S"(A) - 27. (27) 



In the case of the neo-Hookean model Q and thus in the case of the 
third-order elasticity model ([s]) in pure bending, these quantities are 



7 = ixX~ 



a = fiX^ 



2(3 = fi{X^ + A 



-2n 



(28) 



These quantities have already been computed in (Destrade and Scott, 2004; 
Destrade and Ogden, 2005) for the Gent material in general. Specialized to 
the case of plane strain they reduce to: 



a 



fiJmX'^ 



Jm + 2 — A^ 



A 



-2' 



7 



HJmX 



Jm + 2 — A^ — A 



2{Jm + 2 - A2 - A-2) 



A^ + A-2 + 



2(A2-A-2)2 



A 



-2 



For the Fung material in plane strain they are: 



a 



fiX e 



2pfe(A2+A-2-2) 



7 = ^A" 



-2p6(A2+A-2-2) 



(3='^[X' + A-2 + 2biX' - A-2)2]eKA^+^-^-2). 

Now if the incremental equations of motion can be solved, subject to the 
following boundary conditions of traction-free inner and outer faces: 







at 



(29) 
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then an infinity of equilibrium states exists adjacent to tlie large bending, 
signaling the onset of instability (in the linearized sense). When this solution 
is determined, we call Acr = ooTa the value of the azimuthal stretch A2 on the 
inner face. It is the critical stretch of contraction of azimuthal line elements. 
Simple calculations show that 

uA = (A-/ - Xi)/4, (30) 

which allows for the complete determination of the current (deformed) geom- 
etry, just prior to instability. In particular, the angle of bending, and inner 
and outer radii follow as 

if = 2iuA)^, ^ = %, T = ^^' (31) 
^ ^ M A uA A wAAcr 

respectively. Note that if it turns out that (f > 27r, then the conclusion is 
that the block can be completely bent into a cylinder without encountering 
any instability phenomenon. 



4 Numerics 



The direct numerical resolution of the differential system (24) is not easy to 
implement because of numerical stiffness issues, especially for thick blocks. 

First, we use the compound matrix method (Haughton and Orr, 1997) 
which is rendered optimal here by the use of the Stroh formulation. It gives 
direct access to the azimuthal critical stretch of contraction. 

Another option is to integrate numerically the non-linear Riccati equation 
satisfied by the impedance matrix, see (Biryukov, 1985). Here we use that 
integration to compute the full mechanical fields inside the bent block. 

4.1 Critical stretch 

Let rj^^\r) and r}^'^\r) be two linearly independent solutions to the equa- 



tions of equilibrium ( 24 ) . Use them to calculate the following six compound 
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variables 0,, 



(1) (2) 
(1) (2) 

(1) (2) 

mm 

(1) (2) 



(1) (2) 

Vi m 

(1) (2) 
^3 ^3 

(1) (2) 

mm 

(1) (2) 

mm 



(1) (2) 
(1) (2) 

mm 

(1) (2) 

mm 

(1) (2) 

mm 



(32) 



Now compute their derivatives with respect to r and find that they satisfy 
1 



(f){r) = -A{r)(f){r), where 



>1,<P2,(P3,<P4,<P5,' 



(33) 



and A[r) is called the compound matrix. Here we find that its non-zero 
entries are 



A21 = 


A5i = A 


A42 = 


-A23 = 


A32 = 


-A24 = 


Am = 


-A33 = 


At r 


= a, we 



65 — Aq2 — —^12, 

A45 = -^53 = n{l - ai/a), 

A35 = -A54. = n, 
2 — (Ti/a, 



Am 
A41 

Ai3 

A31 



-All = CTi/a, 

^63 = /^ll; 

A4Q = -l/a, 

A&4 = —K22- 



(34) 



traction), so that 02 = 03 = </ 
non-zero: this is the initial value 



■- Sre = (inner face free of incremental 
'4 = 05 = 06 = there, and only 0i is 



0(a) = 0i(a)[l, 0,0, 0,0,0]*. 



(35) 



At r = 6, the tractions must also vanish so that 06 must be zero there: this 
is the target condition, 

Mb) = 0. (36) 



Integrating numerically the initial value problem (33)-(35) poses no dif- 
ficulty. Note that the Stroh formulation yields a most simple and optimal 
form for the elements of the compound matrix (other formulations, for ex- 
ample (Haughton, 1999; Coman and Destrade, 2008) involve derivatives of 
the elastic moduli with respect to r). 

Once the dimensions of the block are given, and its strain energy density 
is fixed, it remains to adjust A so that (36) is satisfied. This search yields the 
critical stretch of contraction for circumferential line elements on the inner 
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face and by extension, the critical angle of bending, see previous section. 
Note that to each mode number p corresponds a different value of the critical 
stretch; however only the highest value is meaningful, as the others cannot 
be reached once the buckling has occurred. 

Take for instance the neo-Hookean model Q. Figure [3] shows the de- 
pendence of the critical stretch on the width-to-thickness ratio L/A. Each 
curve corresponds to a different mode number, but only a small part near its 
maximum is relevant, as highlighted by the thicker stroke. To generate Fig- 
ure [3} the equations were non-dimensionalized, and in the end, only Acr and 
L/{pA) remained as unknowns. Then each curve in the figure is drawn by 
taking p = 1, 2, 3, ... in turn and becomes a scaled version of the p = 1 curve. 
As a consequence, the maximum of each curve is the same. This scaling ex- 
plains the general impression given by the figure: although "dispersion" is 
introduced due to the characteristic dimensions of the block, the value of 
the critical stretch itself is quite insensitive to mode numbers, and remains 
in the neighborhood of 0.562 approximatively (Haughton, 1999; Coman and 
Destrade, 2008). 

For a specific example, take a block with aspect ratio L/A = 3.0. When 
it is modeled by the neo-Hookean, Mooney-Rivlin, or general third-order 
elasticity theory, the numerical calculations of the compound matrix method 
predict that it buckles in bending when A2 = = 0.5613, in mode p = 7, 



see Figure Km a). Then the formulas (31 ) give the bending angle as ip = 246 



and the inner and outer radii as = 0.78A, and rt = 2.5A, respectively, see 
Figure Isl 



4.2 Mechanical fields 

The compound matrix method is most appropriate to obtain quickly and 
accurately the critical stretch of contraction — the "eigenvalue" of the Stroh 



problem (24). By establishing links between the equations satisfied by the 
compound variables and those satisfied by the mechanical field variables, 
Haughton (2008) shows that it is possible to determine those latter fields 
throughout the bent block — the "eigenvectors". Here we present a self- 
contained, alternative protocol for finding the "eigenvalues" and the "eigen- 
vectors", based on the so-called impedance matrix. This approach can be 
dated back to Biryukov (1985), see Fu (2005). 

First, we follow Shuvalov (2003) and call 7Vf(r, r^) the matricant solution 
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Figure 3: Bending instability of rubber: for the neo-Hookean model, the 
Mooney-Rivlin model, and the general incompressible model of third-order 
elasticity, the circumferential line elements on the inner face can be con- 
tracted by 44% at most. Each different plot corresponds to a different mode, 
here the first eight. Each plot can be deduced from the first one by a scaling 
factor. Figure (a) reveals that a block with aspect ratio L/A = 3.0, say, is 
going to buckle in mode p = 7. Figure (b) shows that bent block just prior 
to buckling. 
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to (24); it is defined as the matrix such that 



T7(r) = M{r,ra)'n{ra), M{ra,ra) 
and it has the following block structure 



■(4), 



M(r, To) 



'Mi{r,ra) M2ir,ra)' 



(37) 



(38) 



say. We call S = [Srr,SreY the traction vector, and U = [U,VY the dis- 
placement vector, so that 77 = [U, irSy. We use the incremental boundary 
condition Sivr, 



in (|37|) and (jSSj) to find that 
rS{r) = z°'{r)U{r), where z°- e 



-1 



(39) 



is the conditional impedance matrix (Shuvalov, 2003) (here "conditional" 
refers to the assumed inclusion of the traction-free incremental boundary 
condition at r = Tq). 

Substituting the impedance matrix z"- into the incremental equations of 
equilibrium (24) gives 

d 



dr 



U 



GiU 



Goz'^U, 



dr 



[z'^U) = ^GsU + ^-G+z^U, (40) 



where Gi, G2, G3, and Gf = G^ are the 2x2 sub-blocks of G. Eliminating 
U between these two equations results in the following Riccati differential 
equation for z"'. 



dr 



0, 



(41) 



where the initial condition follows from (39)2 and (37)2- Note that this 
Riccati equation is real because G2 and G3 are Hermitian, see (25), and so 
is z'^ (Shuvalov, 2003). 



Now integrate (|41j) numerically, and adjust the azimuthal stretch ratio so 
that the incremental boundary condition of a traction-free face is satisfied on 
r = Tfe. In other words, find the critical value Acr for A2 so that det z°'{rf)) = 0. 
Then S{rb) = z°'{rh)U{rh) = means that 



(42) 
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This ratio determines the shape of the axial creases on the outer face of the 
bent block. Note however that their amplitude is not known because the 
stability analysis is linear (that is, linearized in the neighborhood of a large 
bending). 

In principle we should be able to integrate simultaneously (40)i and (41 )i 
from Th to r in order to determine the displacement field. In practice, this 
computation can run into numerical difficulties, by encountering singulari- 
ties (Biryukov, 1985). Instead, we use z^, the other conditional impedance 
matrix, found from the condition that S{rb) = 0. We thus integrate simul- 
taneously 

= -GiU--G2z''U, ^z'' = - \z''G2z'' + Gs - iz^'G, + iG+z''] , 

(43) 

with initial conditions 

U{n) = U{n) [1, ~z^nin)/zUrb)f , z\n) = O. (44) 
This procedure is robust, and it gives the entire mechanical displacement field 



throughout the thickness of the bent block. Equation (39)i then gives the 
incremental stress field. Alternatively, the stress field can be found from the 
Riccati differential equations satisfied by the conditional admitance matrices 
(Shuvalov, 2003). 

Figure |4] shows the L/A = 3.0 rubber block of the previous section, in its 
incrementally buckled state. Ten circumferential lines are displayed, clearly 
showing the predicted seven axial wrinkles, and the strong localization of 
the displacement near the inner bent face. In particular, we find that the 
displacement amplitude on that face is more than 3000 times the amplitude 
on the outer face. 



4.3 Results for strain-stiffening models 

In their study on the stability of compressed blocks, Gent and Cho (1999) 
conjecture that it is "not generally necessary to consider stress-strain re- 
lations incorporating finite-extensibility effects", at least for normal Gent 
rubbery materials, for which Jm lies between 20 and 200. Indeed, using 
the numerical techniques exposed in the previous sections, we find that the 
bifurcation curves for the Gent material ([T]) at = 200 are virtually in- 
distinguishable from those of the neo-Hookean material (which corresponds 
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Figure 4: Incremental buckling of a bent rubber block. Here the solid is 
modeled by the neo-Hookean, Mooney-Rivhn, or general third-order elastic- 
ity solid. Its original length-to-thickness ratio was L/A — 3.0. The theory 
predicts that seven creases appear on the inner face, once the critical angle 
of bending is reached. 
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to Jm = oo). Even at Jm = 20.0, the curves are very close to those of the 
neo-Hookean sohd, raising the average critical stretch ratio from 0.562 to 
only 0.564. 

However the situation is different for values of Jm compatible with the 
tensile behavior of biological soft tissues. Hence for = 2.3, a value mea- 
sured for a young human thoracic aorta (Horgan and Saccomandi, 2003), we 
find that the critical stretch ratio is raised by few percents, at around 0.59. 
For Jm = 0.4, a value for a 70 year old human thoracic aorta (Horgan and 
Saccomandi, 2003), it is raised to 0.73 approximatively, showing that the 
earlier the finite-extensibility effects are felt, the more unstable the block is 
in bending. 

We also find that the stiffer Gent materials buckle with less wrinkles than 
the softer ones. Take again the block with length-to-thickness ratio L/A = 
3.0. For Jm = oo, 7 wrinkles appear when cf) = 246°, see previous section; 
for Jm = 20.0, 6 wrinkles appear, at almost the same angle of bending; for 
Jm = 2.3, we have 4 wrinkles, at = 215°; and for Jm = 0.4, we are down 
to 2 wrinkles, at = 113°. 

These results are reported in Figure |5| 

For the Fung material ([3]), the numerical results are remarkably similar. 
Hence when the model is adjusted to account for the behavior of a 21-year- 
old artery, the stiffening parameter is roughly b = 1.0, and the critical stretch 
ratio of contraction is about Acr = 0.62; when it is adjusted for a 70-year-old 
artery, we find that Acr = 0.72, see Figure [6j 

For other, stiffer, soft tissues, the instability occurs even earlier. Hence we 
conducted tensile tests on porcine skin, and estimated that Jm < 0.1 when it 
is modeled with the Gent material, and that 6 > 20 when it is modeled with 
the Fung material, see Figure [7j Then, the corresponding critical stretch of 
contraction is about A^ = 0.80. 

5 Concluding remarks 

The neo-Hookean material, the Mooney-Rivlin material, and consequently, 
the general third-order elasticity incompressible material are not represen- 
tative of strain-stiffening solids (for instance, they do not stiffen in shear). 
For these classes, often used to model rubbers or gels, the critical stretch 
of contraction for circumferential line elements in bending is Acr = 0.56. In 
contrast, the Gent and the Fung material are strain-stiffening materials par 
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Figure 5: How strain stiffening affects bending instabiiity: Wfien tfie soiid 
is described by tfie neo-Hookean model, tfie Mooney-Rivlin model, the gen- 
eral incompressible model of third-order elasticity, or a "rubber" Gent model 
(20 < Jm < 200), the instability occurs when the circumferential line ele- 
ments are contracted by 44% (Acr = 0.56); When it is described by a "young 
artery" Gent model {Jm = 2.3), they can contract by 41% (Acr = 0.59); 
When it is described by an "old artery" Gent model (J^ = 0.4), they can 
contract by 27% (Acr = 0.73) only, before bending creases form. The figures 
on the right show how much a corresponding block with L/A — 3.0 can be 
bent before buckling. 
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Figure 6: How strain-stiffening affects bending instability: Here the solid is 

described by the Fung exponential model for arteries. For a "young artery" 
Fung model (6 = 1.0), the circumferential line elements can contract by 38% 
(Acr = 0.62); For an "old artery" Fung model (6 = 5.5), they can contract by 
^ ('^cr = 0.72) only, before bending creases form. 
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Figure 7: Tensile test on pig skin: nominal stress versus stretch ratio. 
Experimental data and curve fitting using the Gent model (on the left) and 
the Fung model (on the right). Curve fitting analysis gives Jm — 0.08 with 
the Gent model and b ~ 21.5 with the Fung model. 

excellence. When they are used to model stiff (old) arteries, the theory of 
incremental instability gives Acr — 0.73. For pig skin, it is even raised further, 
to about Acr — 0.8. We may thus conclude that the strain-stiffening effect 
actually promotes bending instability. 

Of course, these are theoretical and numerical predictions. Nevertheless, 
it must be kept in mind that the Gent and the Fung materials are popular 
models in the biomechanics literature and in bioengineering simulations. It 
must also be remembered that, irrespective of how advanced a Finite Element 
Analysis software package is, and of how precisely the image and geometry 
of a given soft tissue are captured, numerical simulations are, in fine, as good 
— or as bad — as the constitutive equations they rely upon. 

If a biological soft tissue exhibits strain-stiffening effects, then these must 
be integrated into the constitutive model. If the model predicts instabilities 
when none are observed, then the model must be refined or abandoned. In 
fact different models react differently to different types of instabilities, and 
trends or classifications of models with respect to instabilities are hard to 
decipher in non-linear elasticity (Goriely et al., 2006). 

Another major effect exhibited by biological soft tissues is that of aniso- 
tropy, due to the presence of collagen fiber bundles embedded into the elastin 
matrix. The introduction of preferred directions into a constitutive model can 
have dramatic repercussions with respect to stability analysis (Destrade et 
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al., 2008) and it seems crucial to be able to incorporate and quantify this ef- 
fect. Other important differences between elastomers and soft tissues include 
non-linear viscoelasticity, heterogeneity, growth, and spatial distribution of 
fibers. 
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